Trapped Fermi gases with Rashba spin-orbit coupling in two dimensions 
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We use the Bogoliubov-de Gennes formalism to analyze harmonically trapped Fermi gases with Rashba- 
type spin-orbit coupling in two dimensions. We consider both population-balanced and -imbalanced Fermi 
gases throughout the BCS-BEC evolution, and study the effects of spin-orbit coupling on the spontaneously 
induced countercirculating mass currents and the associated intrinsic angular momentum. In particular, we 
find that even a small spin-orbit coupling destabilizes Fulde-Ferrel-Larkin-Ovchinnikov (FFLO)-type spatially 
modulated superfluid phases as well as the phase-separated states against the polarized superfluid phase. We 
also show that the continuum of quasiparticle and quasihole excitation spectrum can be connected by zero, one 
or two discrete branches of interface modes, depending on the number of interfaces between a topologically 
trivial phase (e.g. locally unpolarized/low-polarized superfluid or spin-polarized normal) and a topologically 
nontrivial one (e.g. locally high-polarized superfluid) that may be present in a trapped system. 

PACS numbers: 05.30.Fk, 03.75.Ss, 03.75.Hh 



I. INTRODUCTION 

The coupling between a quantum particle's intrinsic angu- 
lar momentum (spin) and its center of mass (orbital) motion 
has important consequences in a variety of modern condensed 
matter problems, ranging from quantum spin Hall systems to 
topological insulators and topological superconductors 
This interaction is referred to as the spin-orbit coupling, and 
it arises from coupling of the electron's spin to the local mag- 
netic field that is induced in the electron's reference frame, 
due to the time varying electric field produced by the charged 
background. Since both the strength and the symmetry of the 
spin-orbit coupling are mainly determined by the electronic 
structure of the crystal in condensed matter systems, it is more 
desirable to engineer spin-orbit coupling in alternative sys- 
tems that allow more experimental control over its parame- 
ters. Given the recent experimental advances in simulating 
artificial gauge fields with neutral quantum gases JH-@], it is 
arguable that the prime candidate for engineering spin-orbit 
couplings in a controllable many-body setting seems to be the 
atomic ones. For instance, this has recently been achieved first 
with bosonic yl |4j] and then fermonic ||5l |6[] atomic gases, by 
coupling the momentum of atoms to their spin, via Raman 
dressing of atomic hyperfine states with a pair of laser beams. 
While the symmetry of all of the experimentally engineered 
spin-orbit couplings is so far an equal mixture of Rashba and 
Dresselhaus types, theoretical proposals for creating unequal 
combinations are also underway. 

Since the realization of spin-orbit coupled BECs [3], there 
has been growing theoretical interest in studying spin-orbit 
coupled Fermi gases, even prior to their very recent realiza- 
tion 05|,|6j]. For population-balanced uniform systems, it has 
been shown that the BCS-BEC evolution is a crossover, and 
this evolution can be driven either by increasing the inter- 
particle interaction strength for a fixed spin-orbit coupling or 
by increasing the spin-orbit coupling for a fixed interaction 
strength (no matter how small the interaction strength is) (0- 
[l2tl . On the other hand, for population-imbalanced uniform 
systems, the BCS-BEC evolution is not a crossover, and quan- 
tum phase transitions are found between thermodynamically 



stable and topologically distinct gapped and gapless super- 
fluid phases. These phases are distinguished in momentum 
space by their numbers of zero-energy points, rings or sur- 
faces (depending on the type of spin-orbit coupling) in their 
quasiparticle/quasihole excitation spectrum IU3I419II . 

In direct application to atomic systems, the thermodynamic 
phase diagrams obtained in these works can be easily used 
to extract information about the trapped Fermi gases, at least 
within the semiclassical local-density approximation I20I - I23I1 . 
This commonly used approximation works better and better 
when the number of fermions is increased towards infinity, 
as the finite-size effects become negligible. However, a fully 
quantum mechanical method, e.g. Bogoliubov-de Gennes 
(BdG) formalism, suits better for studying finite-size effects. 
Therefore, in this paper we develop a self-consistent BdG for- 
malism to study harmonically trapped Fermi gases with spin- 
orbit coupling. We only consider the Rashba-type spin-orbit 
coupling in two dimensions due to its numerical simplicity 
(see Sec. Ill Al l, and hope that some of our qualitative conclu- 
sions hold in three dimensions as well. However, we note 
that, given the recent realization of two-dimensional Fermi 
gases ll24l,l25ll . it may also be possible to engineer spin-orbit 
coupling in reduced dimensions. Our main focus here is about 
the spin-orbit coupling induced countercirculating mass cur- 
rents, where we systematically analyze their dependence on 
the spin-orbit coupling, two-body binding energy and popu- 
lation imbalance. We note that induced currents in trapped 
atomic systems have recently been discussed for an optical 
lattice model 12611 . While the Hamiltonian used and the BdG 
formalism developed in this work is completely different, our 
results are in qualitative agreement with each other when there 
is an overlap. 

The rest of the manuscript is organized as follows. In 
Sec. UH we generalize the BdG formalism to spin-orbit cou- 
pled Fermi gases, and derive the self-consistency (order pa- 
rameter and number) equations, probability current density, 
and the associated angular momentum. These equations are 
numerically solved and analyzed in Sec. [TTTl and our main 
findings are briefly summarized in Sec. [IV] 
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n. BOGOLIUBOV-DE GENNES FORMALISM 

Our analysis is based on the self-consistent BdG formal- 
ism, which enables us to include the single-particle quantum 
harmonic oscillator solutions exactly in the real space mean- 
field Hamiltonian density. For this purpose, let us first de- 
scribe the generalization of this theoretical framework to two- 
dimensional trapped Fermi gases with Rashba-type spin-orbit 
coupling. 



A. Hamiltonian and self-consistency equations 

To describe the spin-orbit coupled Fermi gases 
with attractive and short-range interactions, we use 
the Hamiltonian density (in units of h = ks = 1), 
H(r) = E ffl ^J«^W^W + A(r)^j(r)^(r) + 
A*(r)V>4.( r )V't( r )j where the operators i\y a (r) and 
ip<r( r ) create and annihilate a pseudo-spin a fermion 
at position r, respectively, and A(r) is the mean-field 
superfiuid order parameter. The diagonal operator 
K aa {r) = -V 2 /(2M) - p a + V(r) includes both 
the kinetic energy and the harmonic trapping potential 
V(r) = Mw 2 r 2 /2, where M is the mass and [i a is the chem- 
ical potential of a fermions, and uj is the trapping frequency. 
The off-diagonal operator K^{r) = K^Jr) = a(p y + ip x ) 
is the Rashba-type spin-orbit coupling, where a > is 
its strength and pj — ~id/dj is the momentum operator. 
In the polar coordinate system (r, 9), this term becomes 
K^(r) — e~ l6 [d/dr ~ id/(rd0)], which makes the Rashba- 
type spin-orbit coupling numerically much easier to simulate 
in two dimensions due to its rotational invariance. 

The mean-field Hamiltonian can be diagonalized via a gen- 
eralized Bogoliubov-Valatin transformation, and the resultant 
BdG equation can be written as H(r)(p n (r) = e n ip n (r), 
where 
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is the Hamiltonian matrix given in the (p n ( T ) — 
[u t „(r),it in (r),u t „(r),w in (r)] T basis, and s n > are 
the energy eigenvalues. The mean-field superfiuid or- 
der parameter A(r) = g(^(r)ipi(r)) , where g > 
is the strength of the attractive interaction between j" and 
I fermions, and (• • • ) is the thermal average, becomes 
A(r) = gErMn(r)vlJr)f(-e n )+u in (r)v* tn (r)f(e n ). 
Here, f{x) — l/(e x ^ T + 1) is the Fermi function and T 
is the temperature. We may relate g to the energy £& < 
of the two-body bound state between an j" and a I fermion 
in vacuum via the relation, 1/g = 5Z k l/(2£ic — 
where — k 2 /(2M) is the kinetic energy. This leads to 
g — 47r/[Mln(l + 2e c /|Et|)], where e c is the energy cut- 
off used in the k-space integration (e c is specified below in 
Sec. lITll i. The order parameter equation has to be solved self- 
consistently with the number equations N a = j drn a (r), 



where n (T (r) = (r)^ CT (r)) is the local density of a 
fermions. Using the Bogoliubov-Valatin transformations, we 
obtain n a {r) = £„[Kn(r)| 2 /(£ n ) + K„(r)| 2 /(-e„)]. 
Thus, the order parameter and number equations form a closed 
set, determining A(r) and \i a for any given eb, a and T. 

We take advantage of the rotational invariance of 
the Hamiltonian, and conveniently expand the normal- 
ized wave functions as u t „(r) = £) n Cf mn <f> nm (r) 
and Vf n (r) = ^2 n df mn 4> n ,m+i(r) for the t compo- 
nents, and u± n (r) = J2n c imn(f>n, m +i{r) and v± n (r) = 
En rf M^mW for the | ones. Here, <f> nm {r) = 
Rnm{T)Q m {6) are the solutions for the single -particle 
quantum harmonic oscillator problem, where i? ram (Y) = 

/3l" l l+V2«!/(«+ H)!e-^ r2 / 2 rl m 'Lk m| (/3 2 r 2 ) is the ra- 
dial, and Q m {6) — e" n9 /\/2tt is the angular part of the wave 
function. The quantum numbers n = 0, 1, 2, • • • , oo and 
rn = 0, ±1, ±2, • • • , ±oo correspond, respectively, to the ra- 
dial and angular degrees of freedom, /3 = y/Moj, and ill"' (x) 
is the associated Laguerre polynomial. This particular choice 
allow us to decouple the BdG equations into independent sub- 
spaces of m sectors as shown below. 

Using the orthonormality relations 

rdrR nrn (r)R n , m (r) = 5 nn , and d9Q m (6)G m (0) = 
1, where S nn i is the Kronecker delta, this proce- 
dure reduces the BdG equation given in Eq. O to a 
+ 1) x A(n max + 1) matrix eigenvalue problem, 
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(2) 



for each m sector, if we allow < n < n max 
states (n max is specified below in Sec. Hill) . Here, 

-K-am = [w(2n + |m| + 1) — fi a ]S n n' are the single-particle 
terms, S™' = -a rdrR n . m+1 (r){d/dr-m/r)R n , m (r) 

are the spin-orbit coupling terms leading to = 

-a r o °° rdrR 7hm+1 (r)[/3 2 r + (\m\ - m)/r]R n , m (r) + 

2aj3^fri + \m\ + 1 J °° rdrR 7hm+1 (r)R n > ,\ m \ +1 (r), and 
A™' = rdrA(r)R nm (r)R n , m (r) are the pairing terms. 

The same procedure also reduces the order-parameter equa- 
tion to 

A(r) = — [ci mn df mn 'R n m+ i(r)R n i m+1 (r)f(e mn ) 
Zn » — ' 



~t~ C^mnd^rnn' Rnmy) Rn' my) f ( ^-mn)]5 



(3) 



where A (r) = J d?A(r)/(27r) is averaged over the an- 
gular direction (recall the rotational invariance of the sys- 
tem) and it is assumed to be real without loosing general- 
ity, and the angular averaged local-density equations n a (r) = 
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J 27r drn CT (r)/(27r) to 

7l-j-(r) = — y \c^mn c \mn' ^nmW^n' m{ r )f i^mn) 
Z7T * — ' 

mnn' 

~t~ d^-jjijid^riin' Rn,?n-\-l \f}Rn' ,m+l (r)/(-e mn )], (4) 

r H( r ) = 7T~ y] [cimnClmn' Rn,m+l(r)Rn> ,m+l(r) f {£mn) 
Z7T 4 — ' 

mnn' 

"I - d^fYiYidj^mn' Rnmij'^) Rn' my ) J V ^mn)]' (^) 

We recall that the sums are only over the quasiparticle 
states with e mn > 0. Using the orthonormality relations, 
we also obtain the total number of a fermions as N a — 
Emn[4mn/(£mn) + d 2 amn f {-e mn )]. We emphasize that 
these mean-field equations can be used to investigate the low 
temperature properties of the system for all values of Eb and a, 
but they provide only a qualitative description of the system 
outside of the weak-coupling regime, i.e. in the BCS-BEC 
evolution. 



spin-orbit coupling induced particle flow. The angular mo- 
mentum is along the z direction, and its density l„(r) can 
be shown to be related to the strength of the current den- 
sity via £ a (r) — MrJ a (r). Using the orthonormality rela- 
tions, we obtain the total angular momentum of a fermions 
L a = J dr£ a (r) as 

L t = [ mC 1-mnf(£mn) ~ (TO + l)d 2 >m f (-S mn )] , (8) 

mn 

L i = [( m + 1 ) c |mn/( £ mn) ~ md^J {-£ mn )] ■ (9) 
mn 

Having generalized the theoretical BdG framework for the 
trapped two-dimensional Fermi gases with Rashba-type spin- 
orbit coupling, next we discuss our numerical results that 
comes out of this formalism. 



III. NUMERICAL RESULTS 



B. Countercirculating mass currents 

Once the quasiparticle energies and the corresponding wave 
functions are obtained, through self-consistently solving the 
BdG equations discussed above, it is a straightforward task 
to calculate other observables of interest. For instance, next 
we illustrate how we obtain the density of spin-orbit coupling 
induced currents, as well as the intrinsic angular momentum 
associated with the flow of particles. 

Similar to the usual a = treatment, the quantum me- 
chanical probability-current operator for a fermions can be 
identified from the continuity equation. While the presence 
of a spin-orbit coupling leads to additional terms in the to- 
tal particle current operator, these terms do not contribute to 
the current since the expectation value (^il(r)^(r)) = 0. 
Therefore, using the Bogoliubov-Valatin transformation, the 
local current density J CT (r) = [l/(2Mi)](^( r )V^(r) - 
H.c.) circulating around the center of the trapping potential 
becomes J CT (r) = [l/(2Afi)] E»[<„(r)Vu CTn (r)/(e n ) + 
u* n (r)Vu CT „(r) /(—£„) — H.c], where H.c. is the Hermitian 
conjugate. Since Jer(r) circulates along the 8 direction, i.e. 
J CT (r) = J a {r)9, we find 

M r ) = 2nMr i m E C ^mnRnm{r)] 2 f{e m n) 

rn n 

- (m+ l)[Y d 1-^nRn, m +i(r)} 2 f(s mn )}, (6) 

n 

Ji(r) = 2 J Mr Y {( W + X )[zC C imnRn, m +l(r)} 2 f(e m n) 
rn n 

- m[y^^di mn R nm {r)} 2 /(-£ mn )} , (7) 

n 

for the strengths of the current densities. 

In this paper, we are also interested in the intrinsic an- 
gular momentum associated with the spontaneous flow of 



In our numerical calculations, we set a large energy cutoff 
e c £p, and numerically solve the self-consistency Eqs. (f2]l- 
©. Here, e F = k%/(2M) = Mu 2 r 2 F /2 is a characteristic 
Fermi-energy scale, where rp is the Thomas-Fermi radius and 
kp is the Fermi momentum corresponding to the total den- 
sity of fermions at the center of the trap when g = 0, i.e. 
n-)-(O) + nj.(0) = kp/(2ir) at r = 0. We also relate the en- 
ergy cutoff and Fermi energy to the occupation of harmonic 
oscillator levels as e c = oj(N c + 1) and £f = ui(Np + 1), 
respectively, where N c 3> Np. This leads to a total of 
N = (Np + l)(N F + 2) fermions, and therefore, e F w ljVW 
when Np 3> 1. In addition, in order to be consistent with the 
energy cutoff, we choose n max — (N c — \m\)/2 as the maxi- 
mum radial quantum number for a given m, and m max = N c 
as the maximum angular quantum number. In particular, here 
we choose Np = 25 and e c = lep, which corresponds to a 
total of N = 702 fermions and iV c = 181. We checked that 
these values are sufficiently high for the parameter regime of 
our interest, since our results for the order parameter and den- 
sity of fermions agree well (within a few percent) with those 
obtained within the local-density approximation. 

Next, we present our numerical results for population- 
balanced (P = 0) as well as -imbalanced (P ^ 0) Fermi 
gases, where P = (Nf — N±)/N is the population-imbalance 
parameter. 



A. Population-balanced Fermi gases 

In Fig.Q] we set P = and |e&| = 0.2ef, and show A(r), 
n a (r) and J a (r) as a function of r, for a number of a val- 
ues. First of all, since spin-orbit coupling increases the low- 
energy density of states, which is similar to what happens in 
the thermodynamic systems Q2L LXOQ , increasing a monotoni- 
cally increases n a (r) near the center of the trap, and as a result 
of which the Fermi gas shrinks. For instance, when a is in- 
creased from to kp/M, the central n a (r) increases by %25, 
going from k F /(4n) to approximately 5kp/(16n). However, 
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symmetry of the parent Hamiltonian. We note that the direc- 
tions of circulating currents are determined by the chirality of 
the spin-orbit coupling, and the j" and I currents would re- 
verse directions if Kfi(r) = a(p y — ip x ) is used. We see 
that Jj.(r) = — J-f(r) has a nonmonotonic dependence on r: it 
gradually increases from as a function of r making a peak at 
an intermediate distance near the edge of the system, beyond 
which it rapidly decreases to 0. The peak value of J a (r) in- 
creases with increasing a, since a nonzero a is what causes 
counter currents to circulate to begin with. In addition, since 
increasing a shrinks the Fermi gas, the radial location of the 
peak moves inwards towards the trap center. 
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FIG. 1 : (color online) Population-balanced (P = 0) Fermi gas. We 
set the two-body binding energy to |e(,| = 0.2ef and show (a) the 
order parameter A(r) (in units of ef), (b) density n CT (r) [in units of 
fc^/(27r)], and (c) probability current distribution J CT (r) (in units of 
fe|. /M) profiles as a function of radial distance r (in units of tf ), for 
a number of spin-orbit coupling strengths a. 



the corresponding A(r) has a nonmonotonic dependence on 
a. We find that the central A(r) decreases slightly until a crit- 
ical value of a ps 0.5kp/M is reached, beyond which A(r) 
increases with increasing a. The increase in A(r) is again 
mainly a consequence of increased density of states. 

As we discuss below, the presence of a Rashba-type spin- 
orbit coupling spontaneously induces countercirculating mass 
currents. This is clearly seen in Fig. [H C X where the j" and J, 
fermions are rotating around the center of the trap in oppo- 
site directions but with equal speed, due to the time-reversal 



FIG. 2: (color online) We set the two-body binding energy to |et| = 
0.2ef, and show the total angular momentum per a fermion L a /N a 
(in units of h = 1) as a function of spin-orbit coupling strengths 
a (in units of kp/M), for both population-balanced (P = 0) and 
-imbalaced (P = 0.5) Fermi gases. 

In Fig. |2] we show the total angular momentum (per par- 
ticle) associated with the particle flow as a function of a. 
When P — 0, we find that L^/N^ = —L^/N^ monoton- 
ically increases from 0, and we expect it to saturate at 0.5 
when a kp/M. (Since our energy cutoff is not sufficiently 
high compared to the energy associated with the spin-orbit 
coupling when a > 1.2, we could not verify this expectation.) 
We note that the angular momentum of rotating atomic sys- 
tems have so far only been achieved indirectly, by observing 
the shift of the radial quadrupole modes. While this technique 
was initially used for rotating atomic BECs II27[ 12811 . it has 
recently been applied to the rotating fermionic superfluids in 
the BCS-BEC crossover l29ll . We believe a similar technique 
could be used for measuring the intrinsic angular momentum 
of spin-orbit coupled Fermi gases, which may provide an in- 
direct evidence for countercirculating mass currents. 

The origin of spontaneously induced countercirculating 
mass currents can be understood via a direct correspondence 
with the p x +ipj,-superfluids/superconductors [ 26 ] . In these p- 
wave systems, the mass current is associated with the chirality 
of Cooper pairs [30], and this is easily seen by noting that the 
chiral p-wave order parameter Ak oc (x ± iy) ■ k, where k is 
the relative momentum of a Cooper pair, is an eigenfunction of 
the orbital angular momentum with eigenvalue ±ft. This ex- 
plains our findings since it can be shown that the order param- 
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eter of Fermi gases with Rashba-type spin-orbit coupling and 
s-wave contact interactions has chiral p-wave symmetry @|. 
However, unlike the chiral p-wave systems which break time- 
reversal symmetry and belong to the topological class of in- 
teger quantum Hall systems, spin-orbit coupled Fermi gases 
preserve time-reversal symmetry just like quantum spin Hall 
systems, and therefore, they exhibit spontaneously induced 
countercirculating j" and J, mass currents. 



B. Population-unbalanced Fermi gases 

Having presented our numerical results for the population- 
balanced Fermi gases, next we discuss the effects of popula- 
tion imbalance on the system. In Fig. [3] we set P — 0.5 and 
|eb| = 0.2ef, and show A(r), n a (r) and J a (r) as a func- 
tion of r, for a number of a values. When a = 0, we see 
that nf(r) — n\,{r) for r < 0.25ri?, n^(r) > 714.(7") ^ for 
0.25t\f < r < 0.8r F , and rif (r) > nj.(r) = for r > 0.8rF- 
Therefore, the central region corresponds to an unpolarized 
superfluid, and the excess spin-polarized j" fermions are ex- 
pelled towards the edge of the system, i.e. paired j" and J, 
fermions and unpaired normal j" fermions are phase separated, 
with a coexistence region (i.e. a polarized superfluid) in be- 
tween. 

For small a ^ 0, we see that the polarized superfluid re- 
gion rapidly expands towards the central region, and the sys- 
tem mostly consists of a polarized superfluid near the center 
of the trap which is phase separated from a spin-polarized nor- 
mal t fermions residing near the edge. For larger a values, the 
spin-polarized j" gas gives its way to the polarized superfluid, 
and the entire system eventually becomes a polarized super- 
fluid beyond a critical a. This happens around a > O.bkp/M 
when P = 0.5 and |e&| = 0-2ef- We note in passing that 
these findings are consistent with the recent works on thermo- 
dynamic phase diagrams HH-EH, where the phase separated 
state was shown to become gradually unstable against the po- 
larized superfluid phase as a increases from 0. 

In addition, these recent works on thermodynamic systems 
showed that, unlike the a = limit where the unpolarized 
superfluid phase is gapped and polarized superfluid phase is 
gapless, a ^ allows the possibility of having a gapped po- 
larized superfluid phase up to a critical polarization, depend- 
ing on the particular value of a Ill4l4l9ll . Therefore, when 
a ^ 0, in contrast to the topologically trivial unpolarized 
and low-polarized superfluid phases, the polarized superfluid 
phase with sufficiently high polarization becomes topologi- 
cally nontrivial, and has gapless quasiparticle/quasihole exci- 
tations. Note in a trapped system that the topologically non- 
trivial locally high-polarized superfluid phase is sandwiched 
between topologically trivial phases (locally unpolarized/low- 
polarized superfluid and spin-polarized normal) for small a. 

The corresponding A(r) are shown in Fig. 0(b). When 
a = 0, we see that A(r) oscillates with multiple sign changes, 
which is reminiscent of FFLO-type spatially-modulated su- 
perfluid phases. Similar to the P = case, for small a / 0, 
the central A(r) decreases slightly until a critical value of 
a ps Q.lkp/M is reached, beyond which A(r) increases with 
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FIG. 3: (color online) Population-imbalanced (P = 0.5) Fermi gas. 
We set the two-body binding energy to |g(,| = 0.2ef, and show (a) 
the order parameter A (r ) (in units of ef), (b) density n a (r) [in units 
of fef./(27r)], and (c) probability current distribution J a (r) (in units 
of k% /M) profiles as a function of radial distance r (in units of rp), 
for a number of spin-orbit coupling strengths a. 



increasing a. More importantly, the spatial modulations of 
A(r) rapidly disappear with increasing a, and A(r) first be- 
comes finite and then gradually increases near the edge of 
the system. This again indicates that the polarized superfluid 
phase expands towards the edge of the system as a gets larger. 
For larger a values, A(r) gradually increases everywhere, and 
it eventually becomes nearly flat for a substantial region of 
the system, except for a small region around the edge. These 
findings suggest that FFLO-type modulated phases, which are 
known to play a minor role in the thermodynamic phase dia- 
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grams when a = 0, becomes irrelevant for sufficiently large 
a. Therefore, our work provides supporting evidence that the 
recent thermodynamic phase diagrams IIF4T - ll9ll . where FFLO- 
type phases were entirely neglected, are qualitatively accurate 
at least within the mean-field approximation. 

Since population imbalance breaks the time-reversal sym- 
metry when P 0, the j" and I fermions again rotates 
(mostly) in opposite directions with unequal speeds. Simi- 
lar to the P = case, we again see that |Jj-(r)| > 
has a nonmonotonic dependence on r, and the peak value 
of J a (r) increases with increasing a. In Fig. [2] we see 
that \L-f\/N^ > L^/Ni increases from nonmonotonically, 
and we again expect \L a \/N a to be bounded by 0.5 when 
a ^> Uf/M. Having analyzed the A(r), n a (r) and J a {r) 
profiles, and L a , next we analyze the quasiparticle/quasihole 
excitation spectrum of the system. 



C. Inner and outer interface modes 

In Fig. 2J we show e mn as a function of m for population- 
balanced and -imbalanced Fermi gases. First of all, we note 
that the spectrum satisfies e mn = — £— m-l,n, which follows 
from the particle-hole symmetry of the parent Hamiltonian. 
When P = and a — 0, it is well-known that the quasiparti- 
cle and quasihole spectrum are separated with an energy gap 
around m « 0. When P — and a ^ 0, it is expected that 
the spectrum splits into two in m space, creating two identi- 
cal energy gaps located at finite m values. Their locations are 
approximately symmetric around m = 0, and this is clearly 
seen in Fig.@|a). For low P ^ the spectrum is similar. 

When P 7^ is sufficiently high and a is small, we show 
in Fig.@Jb) that the continuum of quasiparticle and quasihole 
spectrum are connected by two discrete branches, i.e. inner 
and outer interface modes [23]. This indicates that there must 
be two phase boundaries (interfaces) between a topologically 
nontrivial superfluid phase and a trivial one. In our case, while 
the inner mode occurs at the interface between the locally un- 
polarized or low-polarized superfluid phase existing near the 
center of the trap and locally high-polarized superfluid phase 
existing at some intermediate region, the outer mode occurs 
at the interface between the locally high-polarized superfluid 
phase and locally spin-polarized normal phase existing near 
the edge of the system. However, the energy separation be- 
tween the inner interface modes becomes larger with increas- 
ing a, which causes this branch to move completely into the 
continuum beyond a critical a value. Therefore, for large 
a, the continuum of quasiparticle and quasihole spectrum are 
connected by a single branch of outer interface modes. This 
is clearly seen in Fig. SJc), and it is a direct consequence 
of the disappearance of the inner phase boundary, which ap- 
proximately happens when a > 0.5kp/M, as discussed in 
SecHHB] 
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FIG. 4: (color online) We set the two-body binding energy to = 
0.2ef, and show the single-particle excitation spectrum e mn (in 
units of ef) as a function of angular quantum number m. Here, the 
population-imbalance parameter P and spin-orbit coupling strength 
a are P = and a = 0Mf/M in (a), P = 0.5 and a = 0Ak F /M 
in (b), and P = 0.5 and a = 1.0k F /M in (c). 



IV. CONCLUSIONS 

To conclude, here we studied harmonically trapped Fermi 
gases with Rashba-type spin-orbit coupling in two dimen- 
sions. We considered both population-balanced and - 
imbalanced Fermi gases throughout the BCS-BEC evolution, 
and paid special attention on the effects of spin-orbit cou- 
pling on the spontaneously induced countercirculating mass 
currents and the associated intrinsic angular momentum. 

One of our main findings is that even a small spin-orbit cou- 
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pling destabilizes FFLO-type spatially modulated superfluid 
phases against the polarized superfluid phase. This suggest 
that FFLO-type modulated phases, which are known to play a 
minor role in the thermodynamic phase diagrams when a = 0, 
becomes irrelevant for sufficiently large a. Therefore, we 
provided supporting evidence that the recent thermodynamic 
phase diagrams lfl4l - fl9ll - where FFLO-type phases were en- 
tirely neglected, are qualitatively accurate at least within the 
mean-field approximation. We also found that the phase sep- 
arated state rapidly becomes unstable against polarized super- 
fluid phase as a increases from 0, which is in good agreement 
with recent works on thermodynamic phase diagrams. In ad- 
dition, we showed for population-imbalanced Fermi gases that 
the continuum of quasiparticle and quasihole excitation spec- 
trum can be connected by zero, one or two discrete branches 
of interface modes depending on the particular value of P 



and a. The number of branches is determined by the number 
of interfaces between a topologically trivial phase (e.g. lo- 
cally unpolarized/low-polarized superfluid or spin-polarized 
normal) and a topologically nontrivial one (e.g. locally high- 
polarized superfluid), that may be present in a trapped system. 
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